a) Field of the Invention
This invention relates to an image processing method for correcting distortions in a ground image acquired from a satellite, and especially to a method for orthocorrecting distortions caused by the relief of a ground surface.
b) Description of the Related Art
Firstly, a description will be made about an image acquired by a satellite (satellite image) and the principle of its correction.
FIG. 6 is a simplified pictorial illustration of acquisition of an image by a satellite, and shows the satellite at numeral 1 and also depicts a ground surface 2 and a ground-surface scanning line 3.
In FIG. 6, the satellite 1 acquires an image of the ground surface 2 by detecting electromagnetic waves from the ground surface 2 with a sensor (not illustrated) mounted on the satellite 1 such that the sensor is directed downwardly. When the sensor is a line sensor, such as CCD, commonly employed as a satellite sensor, the sensor is arranged such that cells (sensor elements) which make up the sensor are arrayed in a row extending at a right angle relative to a flying direction L of the satellite 1. A range on the ground surface 2, which is observed at the same time by the sensor, is the ground-surface scanning line 3. The ground-surface scanning line 3 moves together with the satellite 1 in its flying direction L. As the satellite 1 moves, detection data by the individual cells are inputted in predetermined sampling cycles. The input timing of such detection data is the same for all the cells. The data over the entire ground-surface scanning line 3 are equivalent to a single line of satellite image. Accordingly, data equivalent to a single line of satellite image are acquired by a single input of data (single sampling). As the satellite 1 flies on an orbit, the ground surface 2 is continuously observed at constant swaths along scanning lines 3 as described above.
FIG. 7 is a view showing the principle of a correction of a satellite image obtained as described above. FIG. 7 depicts an observed image 4 and a corrected image 5.
In FIG. 7, the image acquired by the satellite as described above is transmitted as an observed image 4 to the ground. In the observed image 4, configurational distortions have been developed by causes such as fluctuations in the satellite orbit and attitude, distortions in an optical system of the sensor, and terrain on the earth""s surface. To make good use of the image acquired by the satellite, it is hence necessary to correct these distortions in the observed image 4 and to convert it into a corrected image in accordance with a predetermined cartographic projection method. This correction is to perform mapping from a pixel A(p,l) in the observed image 4 into a pixel Axe2x80x2(x,y) in the corresponding corrected image 5. This mapping is called xe2x80x9cmapping Fxe2x80x9d.
The observed image 4 is represented by a p-l coordinate system, in which the arrayed direction of the cells in the sensor on the satellite 1 (in other words, the longitudinal direction of the ground-surface scanning line 3) is set as the p-coordinate axis while the traveling direction L of the satellite 1 is set as the l-coordinate axis. The corrected image 5, on the other hand, is represented by an x-y coordinate system, in which a direction corresponding to the p-coordinate axis of the observed image 4 is set as the x-coordinate axis while a direction corresponding to the l-coordinate axis of the observed image 4 is set as the y-coordinate axis.
FIG. 8 is a flow chart illustrating processing for the correction of the observed image 4.
In the corrected image 5 depicted in FIG. 7, pixels are arrayed and set at equal intervals. In the observed image 4, on the other hand, the positions of pixels corresponding to the individual pixels in the corrected image 5 are not arrayed at equal intervals for such distortions as described above. The correction illustrated in FIG. 8 is to find out, with respect to each pixel position in the corrected image 5, a pixel at a corresponding pixel position in the observed image 4 and to place the pixel at the relevant pixel position in the corrected image 5.
In FIG. 8, the correction consists of a distortion modeling step 100 and a resampling step 101. In the distortion modeling step 100, with respect to each pixel position in the corrected image 5, a corresponding pixel position is determined in the observed image 4. Referring back to FIG. 7, the position of the pixel A(p,l) in the observed image 4, said position corresponding to the position of the pixel Axe2x80x2(x,y) in the corrected image 5, is determined using a distortion model 6 prepared in advance. For this purpose, it is necessary to determine an inverse map Fxe2x88x921 of the map F.
Here, the distortion model 6 is estimated, for example, from fluctuations in the orbit and attitude of the satellite, distortions in the optical system of the sensor and the like at the time of the capture of the observed image 4, or is obtained from the results of an observation or the like by the satellite.
When the position of the pixel A(p,l) in the observed image 4, said position corresponding to the position of the pixel Axe2x80x2(x,y) in the corrected image 5, is determined in the distortion modeling step 100, the routine then advances to the resampling step 101. In this step, the pixel intensity at the position of the pixel A(p,l) is determined and is used as a pixel intensity at the pixel Axe2x80x2(x,y) in the corrected image 5.
FIG. 9 is an illustration of the resampling step 101 in FIG. 8, and describes the resampling step 101 on the basis of the observed image 4 and corrected image 5 depicted in FIG. 7. Incidentally, numeral 7 designates the center positions of pixels (hereinafter called xe2x80x9cpixel positionsxe2x80x9d) in the observed image 4.
In the distortion modeling step 100 described with reference to FIG. 8, the position of the pixel A(p,l) in the observed image 4, said position corresponding to the position of the pixel Axe2x80x2(x,y) in the corrected image 5, was determined. As shown in FIG. 9, however, the thus-determined position of the pixel A(p,l) generally does not register with a pixel position 7 in data observed by a sensor. It is, therefore, necessary to determine the intensity of the data observed at the position of the pixel A(p,l) by interpolation from the pixel intensity or intensities at one or more pixel positions 7 actually obtained around the position of the pixel A(p,l). As an interpolating method for this purpose, an interpolation formula such as nearest neighbor, bilinear or cubic convolution is commonly employed.
As a first conventional correction method of satellite-acquired image data, said method being based on the above-described principle, specifically as an orthocorrection method, there is known, for example, the method described in Demystification of IKONOS by Thierry Toutin and Philip Cheng, Earth Observation Magazine, 9(7), 17-21, 2000. According to this method, the distortion modeling step 100 in FIG. 8 is represented by the following formula (1):                                                         l              =                                                                                          b                      ⁢                                              xe2x80x83                                            ⁢                                              (                                                                              Nth                            ⁢                                                          -                                                        ⁢                            order                            ⁢                                                          xe2x80x83                                                        ⁢                            polynominal                            ⁢                                                          xe2x80x83                                                        ⁢                            in                            ⁢                                                          xe2x80x83                                                        ⁢                            x                                                    ,                                                      y                            ⁢                                                          xe2x80x83                                                        ⁢                            and                            ⁢                                                          xe2x80x83                                                        ⁢                            z                                                                          )                                                                                                                                                        a                      ⁢                                              xe2x80x83                                            ⁢                                              (                                                                              Nth                            ⁢                                                          -                                                        ⁢                            order                            ⁢                                                          xe2x80x83                                                        ⁢                            polynominal                            ⁢                                                          xe2x80x83                                                        ⁢                            in                            ⁢                                                          xe2x80x83                                                        ⁢                            x                                                    ,                                                      y                            ⁢                                                          xe2x80x83                                                        ⁢                            and                            ⁢                                                          xe2x80x83                                                        ⁢                            z                                                                          )                                                                                                                                                                    p              =                                                                                          d                      ⁢                                              xe2x80x83                                            ⁢                                              (                                                                              Nth                            ⁢                                                          -                                                        ⁢                            order                            ⁢                                                          xe2x80x83                                                        ⁢                            polynominal                            ⁢                                                          xe2x80x83                                                        ⁢                            in                            ⁢                                                          xe2x80x83                                                        ⁢                            x                                                    ,                                                      y                            ⁢                                                          xe2x80x83                                                        ⁢                            and                            ⁢                                                          xe2x80x83                                                        ⁢                            z                                                                          )                                                                                                                                                        c                      ⁢                                              xe2x80x83                                            ⁢                                              (                                                                              Nth                            ⁢                                                          -                                                        ⁢                            order                            ⁢                                                          xe2x80x83                                                        ⁢                            polynominal                            ⁢                                                          xe2x80x83                                                        ⁢                            in                            ⁢                                                          xe2x80x83                                                        ⁢                            x                                                    ,                                                      y                            ⁢                                                          xe2x80x83                                                        ⁢                            and                            ⁢                                                          xe2x80x83                                                        ⁢                            z                                                                          )                                                                                                                                                    (        1        )            
This formula is called xe2x80x9cthe rational polynominal methodxe2x80x9d. Whenever an image is acquired, the coefficients in the polynominals in x, y and z are calculated and determined by using ground control points GCP as will be described hereinafter. By the formula (1) making use of the coefficients so determined, coordinate positions in the observed image 4, said pixel positions corresponding to the individual pixel positions in the corrected image 5, are determined.
FIG. 10 is an illustration of the first correction method making use of such ground control points.
In FIG. 10, a plurality of ground control points Gi (i=1, 2, . . . ), for example, three ground control points G1,G2,G3 are chosen in an observed image 4. These ground control points make it possible to clearly show shapes and positions in the observed image 4, and position coordinates in the corrected image 5, said position coordinates corresponding to the positions of these ground control points, can be correctly obtained by reading them from a map or by a survey. In the corrected image 5, the positions corresponding to these ground control points G1,G2,G3 are designated as G1xe2x80x2,G2xe2x80x2,G3xe2x80x2, respectively.
Once the plural ground control points Gi(pi,li) in the observed image 4 and their corresponding positions Gixe2x80x2(xi,yi,zi) in the corrected image 5 are determined as described above, the unknown coefficients in the polynominals represented by the formula (1) can be determined by the least-square method from the pixel coordinates of these ground control points Gi(pi,li) and map coordinates of the positions Gixe2x80x2(xi,yi,zi). As a result, inverse mapping Fxe2x88x921 is feasible. Incidentally, the order of each polynominal should be set as needed depending on the number of the ground control points GCP and the extents of distortions.
This first method has a merit in that, insofar as elevation data DEM (digital elevation model) and ground control points in an observation target region are available, distortion modeling and correction processing making use of the distortion modeling are feasible even if neither the orbit data or attitude data of a satellite nor an optical model of a sensor of the satellite is available.
As a second conventional correction method of satellite-acquired image data, there is also known a so-called non-orthocorrection method, namely, a method that does not take into consideration a relief of elevations on the ground surface. This is a method relying upon a geometric model, which makes use of the orbit data and attitude data of a satellite and an optical model of a sensor of the satellite. This method is characterized in that a corrected image is divided into blocks by lattices of equal intervals and distortion model formulas are determined for the individual blocks, respectively.
FIG. 11 is a view for making an explanation on the division into blocks by the lattices in the second correction method.
In the corrected image 5 of FIG. 11, lattices 8 are set at equal intervals in both vertically and horizontally such that the corrected image 5 is divided into plural blocks 9 by the lattices 8. Pixel positions in the observed image 4 are put into groups by these blocks 9. The pixel positions in the observed image 4, which correspond to individual lattice points 10 where these lattices 8 intersect, are determined by inverse mapping Fxe2x88x921, respectively. Lattices 13 are then set with lattice points 12 located at the thus-obtained pixel positions, respectively, so that the observed image 4 is also divided into blocks 11 by these lattices 13. By way of example, FIG. 11 illustrates that a lattice point A in the observed image 4 has been determined as a lattice point Axe2x80x2 in the corrected image 5 by the inverse mapping F31 1 and that these lattice points Axe2x80x2, A correspond to each other.
Every block 11 in the observed image 4 corresponds to any one of blocks 9 in the corrected image 5. Between the observed image 4 and the corrected image 5, the blocks 11 and the blocks 9 correspond to each other on a one-to-one basis. Incidentally, the blocks 9, 11 are set such that a distortion becomes approximately even in each block. Between adjacent blocks, however, distortions may be different. However, the inverse mapping Fxe2x88x921 is set such that distortions continue at a boundary between each two adjacent blocks.
Here, bilinear functions represented by the following formula (2) are used for the inverse mapping Fxe2x88x921.
l=a0+a1x+a2y+a3xy
p=b0+b1x+b2y+b3xyxe2x80x83xe2x80x83(2)
The formula (2) consists of a 0th-order term in x and y, 1st-order terms in x and y, and a term in xy. Owing to this, distortions are made to continue at the boundary between each two adjacent blocks. If high-order terms of x2 and y2 and higher were included, a discontinuity (step) would occur at the boundary of each two adjacent blocks. To eliminate this discontinuity (step), it is necessary to add many terms of higher orders such that a calculation can be conducted regarding the whole image as a single block. This calculation, however, is very complex and is time-consuming. In this illustrated method, therefore, the images 4,5 are each divided into blocks of such a size as allowing to regard that a distortion is even in each block, and with respect to each block, bilinear functions are determined. By the thus-determined bilinear functions, the inverse mapping Fxe2x88x921 is performed. It is, therefore, possible to perform the inverse mapping Fxe2x88x921 by a simple formula like the formula (2). In this case, the addition of the term in xy as described above is to make distortions continue at the boundary between each two adjacent blocks.
FIG. 12 is a flow chart illustrating a method for calculating the coefficients a0,a1,a2,a3,b0,b1,b2,b3 in the above-described bilinear functions (2).
In FIG. 12, the calculation consists of a lattice-point calculation step 200 and a coefficient calculation step 201.
Describing the calculation on the basis of FIG. 11 as an example, the lattice-point calculation step 200 is to determine a coordinate position A(p,l) in the observed image 4, which corresponds to each lattice point Axe2x80x2(x,y) in the corrected image 5. For this determination, inverse mapping Fxe2x88x921 is needed. As it is generally impossible to directly perform the mapping Fxe2x88x921 without ground control points, the coordinate position A(p,l) is determined using the repetitive convergence computation of geometric calculation for the mapping F.
FIG. 13 is a flow chart illustrating the geometric calculation for the mapping F, and FIG. 14 is a view illustrating an array of some pixel positions in the observed image 4. Further, FIG. 15 is a simplified pictorial illustration for explaining the geometric calculation for the mapping F. The mapping F is to determine a lattice point Axe2x80x2(x,y) which corresponds to each coordinate point A(p,l). In FIG. 13, orbital data 14, attitude data 15, a sensor optical model 16 and pixel position data 17 correspond to the distortion model 6 in FIG. 8.
In FIG. 14, four pixel positions A1(p1,l1),A2(p2,l2),A3(p3,l3),A4(p4,l4) are illustrated. Firstly taking the pixel position A1(p1,l1) as an example, the calculation illustrated in FIG. 13 will be described.
In FIGS. 13 and 15, from the coordinate l1 of the pixel position A1(p1,l1), an acquisition time t of the pixel position A1 is obtained (step 300). From the orbital data 14 of the satellite 1 and the thus-obtained acquisition time t, the position of the satellite 1 (satellite position) in an inertial space at the time of the acquisition of the coordinate l1 is obtained (step 301). From the satellite position and the optical model 16 of the sensor on the satellite 1, the position of the optical center of the sensor (sensor optical center) 19 in the inertial space at the time of the acquisition of the coordinate l1 is obtained (step 302).
Further, from the acquisition time t obtained in step 300 and the attitude data 15 of the satellite 1, the attitude of the satellite 1 (satellite attitude) in the inertial space at the time of the acquisition time t-is obtained (step 303). From the satellite attitude, the optical model 16 of the sensor and the pixel position data 17 (in this computation, the coordinate p1 of the pixel position A1(p1,l1)), a line-of-sight vector 22 to the pixel position A1(p1,l1)xe2x80x94said line-of-sight vector 22 extending through a position 21 of a sensor element (in other words, a cell) in a sensor plane 20, said position 21 corresponding to the coordinate p1xe2x80x94is determined in the inertial space (step 304). Obtained next is an intersection 24 between the line-of-side vector 22 extending from the sensor optical center 19 determined in step 302 and a ground surface model 23 (step 305). This intersection 24 represents the position of the pixel position A1(p1,l1) of the observed image 4 on the ground surface model 23 as indicated in terms of a latitude and a longitude. Based on the intersection 24, a coordinate Axe2x80x2(x,y) in the corrected image 5, said coordinate Axe2x80x2(x,y) corresponding to the pixel position A1(p1,l1), is determined by a calculation making use of predetermined cartographic projection 18. As the ground surface model 23, a rotating ellipsoid is generally employed.
However, the coordinate Axe2x80x2(x,y) determined to correspond to the pixel position A1(p1,l1) as described above does not register with the lattice point 10 set in the corrected image 5 (FIG. 11). Therefore, a pixel position adjacent to the pixel position A1(p1,l1) in FIG. 14, for example, the pixel position A2(p2,l1) is chosen, and with respect to this pixel position, a similar calculation is performed to determine, in the corrected image 5, a coordinate position Axe2x80x2(x,y) corresponding to the pixel position A2(p2,l1). If the thus-determined coordinate Axe2x80x2(x,y) still does not register with the lattice point 10 set in the corrected image 5, a further pixel position is chosen and a similar calculation is also performed. Then by setting a virtual pixel position in between the pixel positions actually observed by the above-selected cells in the sensor, a similar calculation is performed. By changing the location of the virtual pixel position, this calculation is repeated such that the results of the calculations converge at the lattice point 10 in the corrected image 5. The pixel position A(p,l) in the observed image 4, which corresponds to the lattice point Axe2x80x2(x,y) in the corrected image 5, is approximately determined, for example, by four to five calculations or so.
The above-described calculation is performed with respect to every lattice point 10 in the corrected image 5, so that pixel positions A(p,l) in the observed image 4, said pixel positions A(p,l) corresponding to the respective lattice points 10 in the corrected image 5, are obtained. The foregoing is the lattice-point calculation step 200 in FIG. 12.
In the coefficient calculation step 201 in FIG. 12, on the other hand, the coefficients a0,a1,a2,a3,b1,b2,b3 in the bilinear functions represented by the formula (2) are determined for each block 9 in the corrected image 5 shown in FIG. 11. In short, bilinear functions are determined for every block 9.
As shown in FIG. 16, for example, assume that a block 11 in an observed image 4 corresponds to a block 9 in a corrected image 5. It is the coefficient calculation step 201 that determines bilinear functions to perform inverse mapping Fxe2x88x921 of a pixel in the block 11 into the block 9. For this purpose, the coefficients a0,a1,a2,a3,b1,b2,b3 in the formula (2) are determined by using the coordinates of lattice points at the four corners of these blocks 9,11.
Here, the lattice points A1xe2x80x2(x1,y1),A2xe2x80x2(x2,y1),A3xe2x80x2(x2,y2),A4xe2x80x2(x1,y2) at the four corners in the block 9 correspond to the lattice points (pixel positions) A1(p1,l1),A2(p2,l2),A3(p3,l3),A4(p4,l4) at the four corners in the block 11, respectively. Using the coordinate pairs ((x1,y1),(p1,l1)),((x2,y1),(p2,l2)),((X2,y2),(p3,l3)),((x1,y2),(p4,l4)) at the four corners in the blocks 9,11, simultaneous equations represented by the following formula (3), namely,
l1=a0+a1x1+a2y1+a3x1y1
p1=b0+b1x1+b2y1+b3x1y1
l2=a0+a1x2+a2y1+a3x2y1
p2=b0+b1x2+b2y1+b3x2y1
l3=a0+a1x2+a2y2+a3x2y2
p3=b0+b1x2+b2y2+b3x2y2
l4=a0+a1x1+a2y2+a3x1y2
p4=b0+b1x1+b2y2+b3x1y2xe2x80x83xe2x80x83(3)
are solved to obtain the desired coefficients a0,a1,a2,a3,b1,b2,b3. By introducing the thus-obtained coefficients into the formula (2), bilinear functions for these blocks 9,11 are obtained. These bilinear functions express the pixel position in the block 11, said pixel position corresponding to the pixel position in the block 9. If the pixel position in the block 11 as expressed by the bilinear functions does not register with the pixel position actually observed by a cell of the sensor in the block 9, interpolation is obviously performed as described with reference to FIG. 9.
By performing the above-described coefficient calculation step 201 (FIG. 12) with respect to every block 9 in the corrected image 5, bilinear functions for the individual blocks 9 are obtained, thereby making it possible to perform the inverse mapping Fxe2x88x921 from the observed image 4 into the corrected image 5.
This second method has the following merits:
(1) Because different bilinear models are used for the respective blocks, an error propagation caused by the least-squares method does not occur unlike the first method.
(2) Distortion models are continuous at the boundary between each two adjacent blocks.
(3) Owing to the use of the formula (3), the amount of computation in the calculation of position coordinates at the individual points in a corrected image is smaller compared with that required in the case of the polynominals of the formula (1).
(4) Modeling can be performed even if there are no ground control points.
(5) If ground control points are available, the accuracy of the modeling can be improved further by estimating errors, which may be contained in orbital data and attitude data of a satellite, on the basis of the ground control points and correcting the errors.
The first correction method, however, involves the following problems:
(1) Since ground control points are indispensable for modeling, the first correction method cannot be applied where no ground control points are available as in a region where no map or the like is ready for use.
(2) The accuracy of modeling is dependent on the number and arrangement of ground control points. Depending on their number and arrangement, the accuracy of estimation of the coefficients in polynominals may deteriorate, and especially, there is a potential problem that an error may propagate through points other than the ground control points.
(3) A calculation by the formula (1) is needed for each pixel of corrected images, leading to an increase in the amount of calculations for the correction.
The second correction method, on the other hand, has problems in that, because approximation is effected with a rotating ellipsoidal model that the ground surface does not roll, the second correction method cannot be applied directly to an orthocorrection which takes into consideration elevation data DEM of an observation target region and ground-relief-associated distortions caused by the relief of the ground.
With reference to FIG. 17, a description will now be made about ground-relief-associated distortions caused by the relief of the ground.
As is illustrated in FIG. 17, there is actually a rolled ground surface 25 as opposed to the ground surface model 23 of the rotating ellipsoid employed in the second correction method. A discussion will now be made about a cell S2 in the sensor plane 20 upon acquiring the observed image 4 by the satellite 10 as described above. A coordinate obtained by expressing an intersection P1 between the line-of-sight vector 22, which extends through the cell S2 and the sensor optical center 19 determined as described above, and the ground surface model 23 in accordance with a predetermined cartographic projection method is set as a pixel position in the corrected image. To tell the truth, however, a coordinate obtained by expressing an intersection P2 between a vertical line, which extends through an intersection Q2 between the line-of-sight vector 22 and the rolling ground surface 25, in accordance with the predetermined cartographic projection method should be set as a correct pixel position in the corrected image. A positional offset on the corrected image, said positional offset being equivalent to the positional offset between P1 and P2, is called a xe2x80x9cGround-relief-associatedxe2x80x9d distortion xcex4xe2x80x9d.
An object of the present invention is to solve such problems, and is to provide a method for orthocorrecting a satellite image without needing ground control points, without causing propagation of an error throughout the image and without making the amount of calculations excessively large while making it possible to correct ground-relief-associated distortions caused by the relief of the ground surface.
To achieve the above-described object, the present invention provides a method for orthocorrecting a satellite-acquired image, including distortion modeling for making positions of pixels in a corrected image and positions of pixels in an observed image correspond to each other and also resampling for interpolating intensities of the pixels in the observed image in accordance with results of the distortion modeling, wherein the distortion modeling comprises the following steps:
(a) setting a plurality of control planes of different elevations relative to a control plane of a zero elevation of a reference rotating ellipsoidal earth model;
(b) dividing the corrected image into plural blocks by lattices of equal intervals;
(c) determining, with respect to each of the pixels in the corrected image, a corresponding one of the blocks;
(d) determining two of the control planes, said two control planes sandwiching an elevation value of the position of the pixel above and below the elevation value in the thus-determined block;
(e) calculating two pixel positions in the observed image, which correspond to the elevation value of the pixel position in the block, by using respective distortion model formulas determined as bilinear functions with respect to the two control planes;
(f) linearly interpolating the thus-calculated two pixel positions by the elevation value of the pixel position in the corrected image to determine the pixel position in the observed image, which corresponds to the pixel position in the corrected image; and
(g) repeating the steps (c) to (f) until pixel positions in the observed image are all determined corresponding to the individual pixel positions in the corrected image.
The method according to the present invention can bring about the following advantageous effects:
(1) Using the given elevation data DEM, ground-relief-associated distortions caused by the relief of a ground surface can be corrected.
(2) Because the different bilinear models are used for the respective blocks, an error propagation caused by the determination of model formulas by the least-squares method does not occur.
(3) Distortion models are continuous at the boundary between each two adjacent blocks.
(4) Modeling can be performed even if no ground control points are available.
(5) If there are ground control points, the accuracy of the modeling can be improved further by estimating errors, which may be contained in orbital data and attitude data of the satellite, on the basis of the ground control points and correcting the errors.
Owing to these advantageous effects, the method according to the present invention for the orthocorrection of a satellite image is practical in both accuracy and the amount of calculations.